Merger of black hole and neutron star in general relativity: Tidal disruption, torus 

mass, and gravitational waves 
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We systematically perform the merger simulation of black hole-neutron star (BH-NS) binaries 
in full general relativity, focusing on the case that the NS is tidally disrupted. We prepare BH-NS 
binaries in a quasicircular orbit as the initial condition in which the BH is modeled by a nonspinning 
moving puncture. For modeling the NS, we adopt the F-law equation of state with F = 2 and the 
irrotational velocity field. We change the BH mass in the range Mbh ~ 3.3-4.6M0, while the 
rest mass of the NS is fixed to be M* = 1.4Mq (i.e., the NS mass Mns « 1.3Mq). The radius 
of the corresponding spherical NS is set in the range _Rns ~ 12-15 km (i.e., the compactness 
GMns/^nsc^ ~ 0.13-0.16). We find for all the chosen initial conditions that the NS is tidally 
disrupted near the innermost stable circular orbit. For the model of TiNS = 12 km, more than 97 % 
of the rest mass is quickly swallowed into the BH and the resultant torus mass surrounding the BH 
is less than 0.04Mq. For the model of iiNS ~ 14.7 km, by contrast, the torus mass is about O.I6M0 
for the BH mass ~ AMq. The thermal energy of the material in the torus increases by the shock 
heating occurred in the collision between the spiral arms, resulting in the temperature 10^*^-10^^ 
K. Our results indicate that the merger between a low-mass BH and its companion NS may form 
a central engine of short gamma-ray bursts (SGRBs) of the total energy of order 10*^ ergs if the 
compactness of the NS is fairly small < 0.145. However, for the canonical values Mns = 1.35M0 and 
Rns = 12 km, the merger results in small torus mass, and hence, it can be a candidate only for the 
low-energy SGRBs of total energy of order 10** ergs. We also present gravitational waveforms during 
the inspiral, tidal disruption of the NS, and subsequent evolution of the disrupted material. We 
find that the amplitude of gravitational waves quickly decreases after the onset of tidal disruption. 
Although the quasinormal mode is excited, its gravitational wave amplitude is much smaller than 
that of the late inspiral phase. This reflects in the fact that the spectrum amplitude sharply falls 
above a cut-off frequency which is determined by the tidal disruption process. We also find that the 
cut-off frequency is 1.25-1.4 times larger than the frequency predicted by the study for the sequence 
of the quasicircular orbits and this factor of the deviation depends on the compactness of the NS. 



PACS numbers: 04.25.Dm, 04.30.-w, 04.40.Dg 



o 



X 



I. INTRODUCTION 

The merger of black hole (BH)-neutron star (NS) bina- 
ries is one of the most promising sources of kilo- meter size 
laserinterferometric gravitational wave detectors such as 
LIGO and VIRGO. Although such system has not been 
observed yet in our Galaxy in contrast to NS-NS bina- 
ries statistical studies based on the stellar evolution 
synthesis suggest that the merger will happen 1-10% as 
frequently as the merger of the NS-NS binaries in the 
universe Because the BH mass should be more than 
twice as large as the canonical NS mass ~ 1.35M0, the 
typical amplitude of gravitational waves from the BH-NS 
binaries will be larger than that from the NS-NS binaries 
even if the distance is larger. This indicates that the de- 
tection rate of the BH-NS binaries by the gravitational 
wave detectors may be more than 10% as high as that 
for the NS-NS binaries, suggesting that the detection of 
such system will be achieved in the near future. 

The final fate of the BH-NS binaries is divided into 
two cases depending primarily on the BH mass: (1) when 
the BH mass is small enough, the NS will be tidally dis- 
rupted before it is swallowed by the BH; (2) when the 



BH mass is large enough, the NS will be swallowed into 
the BH without tidal disruption. The tidal disruption oc- 
curs when the tidal force of the BH is stronger than the 
self- gravity of the NS. Such condition is approximately 
written as 
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where Mns is the NS mass, i?NS is the circumferential 
radius of the NS, Mbh is the BH mass, r is the orbital 
separation, and C is a nondimcnsional coefficient. Using 
the orbital angular velocity Q, Eq. ([T]) may be written as 
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where G is the gravitational constant and Mq = Mbh + 
A/ns. Actually, the latest high-precision numerical study 
for the quasicircular states of the BH-NS binaries shows 
that the tidal disruption occurs if the following condition 
is satisfied [3, Q : 
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where the value of C is « 0.270 for F = 2 poly tropic 
equation of state (EOS). We note that in [1, [1] the NS 
is assumed to be irrotational and the BH spin is set to 
be zero. According to an approximate general relativistic 
study Q , the value of C depends weakly on the stiffness 
of the EOSs for 2 < r < 3, and hence, it would be close 
to 0.27 for the NS which likely has such stiff EOS. 

Substituting the typical values considered in this pa- 
per, the tidal disruption is expected to set in at 
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where c is the speed of the light and q denotes the mass 
ratio q = M^s/Mbh- The corresponding frequency of 
gravitational waves is calculated from / = JI/tt as 
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According to the third post-Newtonian (PN) study, the 
innermost stable circular orbit (ISCO) for q w 1/3 is 
located at GMofl/c^ w 0.11 6]. The tidal effect reduces 
this value shghtly to be GMqVL/c^ - 0.08-0.09 0, i]. 
Adopting this value, we expect that the tidal disruption 
of the NSs of mass 1.3M0 occurs for g > 1/4 if i?NS — 13 
km, and for q ^ 1/3 if i?NS = 11 km. This indicates 
that the tidal disruption occurs only for the binaries of 
low-mass BH and NS at an orbit very close to the ISCO. 

The tidal disruption of NSs by a BH has been stud- 
ied with great interest because of the following reasons: 
(i) Gravitational waves during the tidal disruption may 
bring the information about the NS radius because the 
orbital frequency at the onset of tidal disruption depends 
strongly on the compactness of the NS (GA/ns/^nsc^) 
01 0j 111- The NS mass will be determined by the data 
analysis for observed gravitational waves in the inspiral 
phase . If the NS radius could be determined from ob- 
served gravitational waves emitted during the tidal dis- 
ruption, the resultant relation between the mass and the 
radius of the NSs may constrain the EOSs of the high 
density nuclear matter [1, Q . (ii) The tidally disrupted 
NSs may form a disk or torus of mass larger than O.OIMq 
around the BH if the tidal disruption occurs outside the 
ISCO. Systems consisting of a rotating BH surrounded 
by a massive, hot torus have been proposed as one of the 
likely sources for the central engine of gamma-ray bursts 
with a short duration [H (hereafter SGRBs). Hence, the 
merger of a low-mass BH and its companion NS can be a 
candidate central engine. According to the observational 
results by the Swift and HETE-2 satellites [il|, the total 
energy of the SGRBs is larger than ^ 10^^ ergs, and typ- 
ically lO'^^-lO^" ergs. The question is whether or not the 



mass and thermal energy of the torus are large enough 
for driving the SGRBs of such huge total energy. 

The tidal disruption of NSs by a BH has been investi- 
gated in the Newtonian [T^ and approximately general 
relativistic [H, [l3| simulations. However, the criterion 
of the tidal disruption and the evolution of the tidally 
disrupted NS material would depend strongly on general 
relativistic effects around the BH, and hence, a simula- 
tion in full general relativity is obviously required. In 
the previous paper p^ . we presented our first numer- 
ical results for fully general relativistic simulation per- 
formed by our new code which had been improved from 
the code used for the NS-NS merger [ll,[l3|: We handle 
an orbiting BH employing the moving puncture method, 
which has been recently developed (e.g., [H, [3). We 
prepared a quasicircular state composed of a nonspin- 
ning moving puncture BH and a corotating NS as the 
initial condition. We computed the BH-NS quasicircu- 
lar state by employing our new formalism based on the 
moving puncture framework. From the merger simula- 
tion, we found that a torus with mass ^ Q.IMq is an 
outcome for Mns = 1.3Mo and i?Ns — 13-14 km, and 
for A/bh = 3-4AfQ. 

In the previous work, we assumed that the NSs are 
corotating around the center of mass of the system for 
simplicity. However, this is not very realistic velocity 
field for the NS in close compact binaries [13]. Recently, 
we have developed a new code for computing quasicir- 
cular states of a BH and an irrotational NS. For obtain- 
ing the irrotational velocity field of the NS, one has to 
solve the equation for the velocity potential [2l|. We 
employ the same method as that for computing the NS- 
NS binaries (e.g., [13, H^), which has been used in our 
project for the NS-NS merger [H, [l3|- (For solving the 
basic equations for the initial data, we use the spectral 
method library LORENE developed by the Meudon rel- 
ativity group 24].) In this paper, we perform the sim- 
ulations employing these new quasicircular states as the 
initial conditions. In addition, we modify our simulation 
code for the Einstein evolution equation to the fourth- 
order scheme for more accurate computation. Changing 
the BH mass and NS radius, we systematically investigate 
the dependence of the torus mass formed after the merger 
on A/bh and i?NS • We also compute gravitational wave- 
forms during the tidal disruption and their spectrum. 

The paper is organized as follows. Section II summa- 
rizes the initial conditions chosen in this paper. Section 
HI briefly describes the formulation and numerical meth- 
ods for the simulation. Section IV presents the numerical 
results of the simulation for the merger of BH-NS bina- 
ries. Section V is devoted to a summary. In the following, 
the geometrical units of c = G = 1 are used. 



II. INITIAL CONDITION 

We compute the initial condition for the numerical sim- 
ulation by employing a formulation in the moving punc- 
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TABLE I: Parameters for the quasicircular states of BH-NS binaries. The mass parameter of the puncture, the BH mass, the 
rest mass of the NS, the mass, radius, and normaUzed mass of the NS in isolation, the ADM mass of the system, the total 
angular momentum in units of A/^, the orbital period in units of Mo — A/bh + A/ns, and the compactness of the system defined 
by (Afo^^)^/^ The BH mass is computed from the area of the apparent horizon A as (A/167r)^/^ MoQ, is ^ 0.040-0.041 for all 
the models. 



Model Mp{Mq) 


AfBH(Ar0) M,{Mq) AfNs(M"0) 


q 


7?NS (km) Af./Ki/2 


M(Mq) 


J/Af2 


Po/Mo 


(Afo!^)^''^ 


A 3.906 


3.975 


1.400 


1.302 


0.327 


13.2 


0.150 


5.227 


0.662 


154 


0.118 


B 3.881 


3.950 


1.400 


1.294 


0.327 


12.0 


0.160 


5.193 


0.660 


154 


0.119 


C 3.930 


4.001 


1.400 


1.310 


0.328 


14.7 


0.140 


5.260 


0.663 


152 


0.119 


D 3.255 


3.321 


1.400 


1.302 


0.392 


13.2 


0.150 


4.576 


0.725 


157 


0.117 


E 3.234 


3.300 


1.400 


1.294 


0.392 


12.0 


0.160 


4.546 


0.721 


154 


0.118 


F 4.557 


4.627 


1.400 


1.302 


0.281 


13.2 


0.150 


5.878 


0.612 


158 


0.117 



ture framework which is described in our previous papers 
Tsl. This formulation is slightly different from that in 
S^U, [2^, [2^ , although both formulations are based on 
the confmormal flatness formulation for the three-metric, 
and solve the equations of the maximal slicing condition, 
and Hamiltonian and momentum constraints. (See p7j 
for the simulations preformed employing the initial data 
presented in Q). Although the basic equations adopted 
in this paper are the same as those described in [15l |. 
we slightly change the method for defining the quasicir- 
cular state as follows: In the previous paper, we deter- 
mined the center of mass of the system imposing that the 
dipole moment of the conformal factor at spatial infinity 
is zero. In this paper, the center of mass is determined 
from the condition that the azimuthal component of the 
shift vector (3"^ at the location of the puncture is equal to 
where fl is the orbital angular velocity. Namely, we 
impose the corotating gauge condition at the location of 
the puncture. The reason for this change is that with this 
new method, the curve of the angular momentum J as 
a function of the orbital angular velocity along quasi- 
circular sequences is closer to the curve derived from the 
third PN equation of motion . Furthermore, the rela- 
tion among the orbital angular velocity, ADM mass M, 
and angular momentum agrees with the results obtained 
by a different method 0, Q fairly well; for the mass ra- 
tio and NS radii adopted in this paper, the disagreement 
of the angular momentum for a given angular velocity is 
- 1%. 

As the velocity field of the NS, we assume the irrota- 
tional one because it is believed to be the realistic ve- 
locity field H^l- The density profile and velocity field 
are determined by solving the hydrostatic equation (the 
first integral of the Euler equation] and an elliptic-type 
equation for the velocity potential |2l| . 

The NSs are modeled by a polytropic EOS 

P^Kp^, (6) 

where P is the pressure, p the rest-mass density, k the 
polytropic constant, and T the adiabatic index for which 
we set r = 2. k is a free parameter and in the follow- 
ing we choose it so as to fix the rest mass of the NSs 
to be Af, = 1.4Afo because this value is close to the 



canonical value [l|. We note that the mass, radius, and 
time can be rescaled arbitrarily by changing k: i.e., if 
we change the value of k from ki to K2, these quantities 
are rescaled systematically by a factor of (^2/^1)^/^ for 
r — 2. By contrast, the nondimensional parameters such 
as M]sis/Rns, Mfl, Mk~^/^, and iiNS's""^^^ are invari- 
ant. In this paper. We present the results for the specific 
values of k for the help to the general readers who are 
not familiar with the nondimensional units. 

Because we choose Af* = IAMq, the NS mass is in the 
narrow range A^ns ~ 1.29-1.31Af0 irrespective of the NS 
radius i?NS- The NS radius is chosen to be 12.0, 13.2, and 
14.7 km. These values agree approximately with those 
predicted by the nuclear EOSs in which i?NS ~ 10-15 km 
for AfNS = 1.3Af0 [H]. The value Rm = 12 km agrees 
approximately with that predicted by the stiff nuclear 
EOSs such as the Akmal-Pandharipande-Ravenhall EOS 

m- 

The BH mass is chosen in the range Af bh ~ 3.3- 
4.6Af0. As we described in Sec. I, the tidal disruption 
of the NS occurs for such low-mass BH. 

In the present work, we prepare the quasicircular states 
with Mofl w 0.04. With such initial condition, the BH- 
NS binary experiences about one and half orbits before 
the onset of merger. In order to reduce the spurious ellip- 
ticity and to take into account the nonzero radial velocity 
associated with the gravitational radiation reaction, we 
may need to perform a simulation for several orbits before 
the onset of merger, starting from a more distant orbit 
[13, mi, mi ■ In the present paper, however, we focus pri- 
marily on the qualitative study of the tidal disruption 
events. For such purpose, the adopted initial conditions 
are acceptable. The simulation of the longterm evolution 
for the inspiral phase is left for the future. 

In Table I, the parameters of the quasicircular states 
adopted in this paper are listed. Specifically, we prepare 
six models for clarifying the dependence of the tidal dis- 
ruption events on the NS radius and mass ratio. The 
NSs for models A, B, and C and for models D and E, 
respectively, have approximately the same mass ratio 
q = Mj^s/Mbu, although the NS radii are different each 
other. Comparing the numerical results among these 
models, the dependence of the properties of the tidal dis- 
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ruption event on the NS radius is clarified. The NS radius 
for models A, D, and F and for models B and E are the 
same, although the mass ratios are different each other. 
Thus, comparison of the numerical results for these mod- 
els clarifies the dependence of the properties of the tidal 
disruption event on the mass ratio. 



III. NUMERICAL METHODS 

The numerical code for the hydrodynamics is the same 
as that in where we use a high-resolution central 
scheme with the third-order piece- wise parabolic interpo- 
lation and with a steep min-mod limiter in which the lim- 
iter parameter b is set to be 2.5 (see appendix A of [33|). 
We have already used this numerical code for the simula- 
tion of the NS-NS binary merger [lilll. We adopt the 
F-law EOS in the simulation as 

P = (r - l)pe, (7) 

where e is the specific internal energy and F = 2. 

For solving the Einstein evolution equation, we use the 
original version of the BSSN formalism [3a|. Namely, 
we evolve the conformal part of the three-metric = 
(ln7)/12, the trace part of the extrinsic curvature K, 
the conformal three-metric 7^ = 7~^^'^7ij, the tracefree 
extrinsic curvature Aij = ^^^^^{Kij — K^ij/3), and a 
three auxiliary variable Fi = S^'^djjik- Here, 7^ is the 
three-metric, Kij the extrinsic curvature, 7 = det(7y), 
and K = Kijj''^ . In the previous paper [l^, we solved 
the equation for 7"^/^. Since then, we have learned that 
we do not have to evolve the inverse of the conformal 
factor even in the moving puncture framework because 
we use the cell-centered grid in which the puncture is 
never located on the grid points. 

For the condition of the lapse function a and the shift 
vector we adopt a dynamical gauge condition in the 
following form: 

{dt- f3'd^)\na^ ~2K, (8) 
dtp' = 0.75f\F,+AtdtFj). (9) 

Here, At denotes the time step in the simulation and the 
second term in the right-hand side of Eq. ^ is intro- 
duced for the stabilization of the numerical computation. 

The numerical code for solving the Einstein equation is 
slightly modified from the previous version as follows: (i) 
All the spatial finite differencing including the advcction 
term such as P^dicf) are evaluated with the fourth-order 
finite differencing scheme. For the advection term, an up- 
wind method is adopted as proposed, e.g., in (ii) For 
the time-integration, a third-order Runge-Kutta scheme 
is employed. With this scheme, the stable evolution is 
feasible. Furthermore, the numerical accuracy is signigi- 
cantly improved, and as a result, the convergent results 
can be obtained by a relatively wider grid spacing than 
by that adopted in the previous papers p^ . 



The location and properties of the BH such as the area 
and circumferential radii are determined using an appar- 
ent horizon finder, for which our method is described in 

For extracting gravitational waves from the geometri- 
cal variables, the gauge- invariant Moncrief variables Rim 
in the flat spacetime 13611 has been computed in our se- 
ries of papers (e.g., [l6l [TtI. ^^). From Rim, the plus 
and cross modes of gravitational waves, /i+ and h^, are 
obtained by a simple algebraic calculation. The detailed 
equations are described in [s^] to which the reader 
may refer. In this method, we split the metric in the 
wave zone into the flat background and linear pertur- 
bation. Then, the linear part is decomposed using the 
tensor spherical harmonics, and the gauge-invariant vari- 
ables are constructed for each mode of the eigen values 
{l,m). The gauge-invariant variables of Z > 2 can be 
regarded as gravitational waves in the wave zone, and 
hence, we derive such modes in the numerical simulation. 

In the present work, we also compute the outgoing part 
of the Newman-Penrose quantity (the so-called ■04; see, 
e.g., [1^ for definition of ■04). By twice time inte- 
gration of 2'04 (and with an appropriate choice of the in- 
tegration constants), one can compute the gravitational 
waveforms. We compare these gravitational waveforms 
with those computed from the gauge-invariant wave- 
extraction method. It is found that the wave phases of 
both waveforms agree well throughout the inspiral phase 
to the tidal disruption phase. However, the amplitude 
does not agree well in the early inspiral state; the mag- 
nitude of the disagreement is ~ 20-30%, although the 
disagreement is much smaller during the tidal disruption. 
The reason is that we extract gravitational waves near the 
outer boundary which is located in the local wave zone 
in the present work, and the amplitude is not the correct 
asymptotic amplitude. To obtain the asymptotic ampli- 
tude, it is necessary to perform the extrapolation using 
the data of different extraction radii. One point worthy 
to note is that the convergence of the amplitude of £'■04 
with changing the extraction radius D is much faster than 
that for DRim- This indicates that -04 is a better tool 
than Rim for evaluating the asymptotic waveforms with 
smaller systematic error, in the simulation that the outer 
boundary is located in the local wave zone. Hence, we 
decide that in this paper, we evaluate the waveform, en- 
ergy, angular momentum, and linear momentum fluxes 
analyzing ?/'4. 

We compute the modes of 2 < Z < 4 for ^04, and found 
that the mode of (^, \m\) — (2,2) is dominant, but I — 
|m| = 3 and Z = |to| = 4 modes also contribute to the 
energy and angular momentum dissipation by more than 
1% for the merger of the chosen initial data. 

We also estimate the kick velocity from the linear mo- 
mentum flux of gravitational waves. The linear momen- 
tum flux dPi/dt is computed from the same method as 
that given in [H, [s^. Specifically, the coupling terms 
between I = \m\ = 2 and I = \m\ = 3 modes and between 
I = \m\ = 2 and {I, \m\) = (2,1) modes primarily con- 
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TABLE II: Parameters for the grid in the numerical simula- 
tion. The grid number for covering one positive direction A^, 
the grid number for the inner uniform-grid domain A'o, the 
parameters for the non- uniform grid domain (Ai,^), the grid 
spacing in the uniform domain in units of Mp, the grid num- 
ber covered for the major diameter of the NS (Lns), and the 
ratio of the location of the outer boundary along each axis to 
the gravitational wavelength at t = 0. 



Model 


A^ 


A^o 


Ai 




Ax/Mp L 


Ns/Aa; 


L/A 


A 


225 


150 


30 


20.5 


1/15 


55 


0.86 


B 


225 


150 


30 


20.5 


1/15 


49 


0.86 


C 


225 


151 


30 


20.5 


1/15 


63 


0.86 


D 


225 


155 


30 


20.5 


1/15 


67 


0.74 


E 


225 


150 


30 


20.5 


1/15 


59 


0.82 


F 


225 


155 


30 


22.0 


1/16 


50 


0.80 


Al 


211 


136 


30 


18.35 


1/13.5 


49 


0.86 


A2 


195 


121 


30 


16.4 


1/12 


44 


0.86 


C2 


195 


121 


30 


16.2 


1/12 


50 


0.86 



tribute the linear momentum flux. From the total linear 
momentum dissipated by gravitational waves, 

AP. = /f<,, (10, 

the kick velocity is defined by APi/M where M is the 
initial ADM mass of the system. 



IV. NUMERICAL RESULTS 

A. Choosing the grid points 

In the simulation, the cell-centered Cartesian, {x, y, z), 
grid is adopted to avoid the situation that the location 
of the puncture (which always stays in the z — plane) 
coincides with one of the grid points. The plane symme- 
try is assumed with respect to the equatorial plane. The 
computational domain of —L < x < L, —L < y < L, and 
< z < L is covered by the grid size of (2iV, 2iV, N) for 
x-y-z where L and N are constants. Following [isl. Isst. 
we adopt a nonuniform grid as follows; an inner domain 
of {2No,2No, Nq) grid zone is covered with a uniform 
grid of the spacing Ax and outside this inner domain, 
the grid spacing is increased according to the relation 
of ^tanh[(j — No)/Ai]Ax where i denotes the i-th grid 
point in each positive direction, and Nq, Ai, and are 
constants. Then, the location of the i-th grid, x^{i), in 
each direction is 

{(i-l/2)Ax 1 < j < A^o 

{i^l/2)Ax + iAiAx (11) 
X log[cosh{(i - A"o)/Ai}] i > Ao 

and x^{—i) — —x^{i) i — 1,2, - ■ ■ ,N for x^ = x and y. 



For all the simulations, we chose A^ = 225 and Ai = 30. 
For most of the simulation, we determined the grid spac- 
ing for the inner uniform domain according to the rule 
Ax/Mp = 1/15 where Mp is the mass parameter of the 
puncture. We have found that with this choice of the grid 
size, a convergent result for the BH orbit is obtained. Ac- 
tually, the simulations for the BH-BH binary in the mov- 
ing puncture framework (e.g., [H,[ll,[3H) show that this 
choice is appropriate for obtaining a convergent result in 
the fourth-order scheme. With increasing the mass ratio 
Mp/M^S, however, the grid number of covering the NS 
becomes small for the fixed value of Ax. We require that 
the major diameter of the NS, L^s, should be covered 
by ^ 50 grid points, because a result with small error 
is obtained with such setting (see Appendix A). By this 
reason, for model F, we choose Ax/Mp = 1/16. 

In Table II, we list the parameters for the grid coordi- 
nates. A denotes the wavelength of gravitational waves 
at t = (A = 7r/f2). Aq is chosen so as to put the NS 
in the inner uniform-grid domain initially. ^ is deter- 
mined from the conditions (i) L/X ^ 0.75 and (ii) the 
grid spacing near the outer boundaries is smaller than 
1.05M. Note that A decreases with the time because the 
orbital radius decreases due to the gravitational radia- 
tion reaction. Hence, L is smaller than A at the onset of 
tidal disruption. The gravitational wavelength after the 
merger sets in is expected to be larger than IIM which is 
approximately equal to the wavelength of the quasinor- 
mal mode (QNM) of the BH excited in the final phase 
of the merger. Because the wavelength is covered by at 
least 10 grid points, we may expect that gravitational 
waves are extracted with a good accuracy near the outer 
boundary. 

To see the dependence of the numerical results on the 
grid resolution, we performed test simulations for Ax = 
Mp/12 (models A2 and C2), and for Ax = Mp/13.5 
(model Al). For models A, Al, and A2 and for mod- 
els C and C2, the locations of the outer boundaries are 
approximately the same, respectively (see Table II). We 
found that the numerical results (see Appendix A) show 
an acceptable convergence for the case that the grid spac- 
ing satisfies the conditions (i) Ax > Mp/12 and (ii) that 
the diameter of the NS is covered by more than w 50 grid 
points. As we discuss in Appendix A, the errors for the 
estimate of the mass surrounding the BH and area of the 
apparent horizon at the end of the simulation are within 
- 10%~30% (- O.OI-O.O2M0) and ~ 0.3%-0.5%, respec- 
tively, for such setting. The total energy and angular 
momentum emitted by gravitational waves are also com- 
puted within ^ 10%-20% error. In the previous papers 
[15|, we performed simulations in the second-order accu- 
rate finite-differencing for solving the Einstein equation. 
In such case, the convergence of the numerical results 
with improving the grid resolution is not very fast. In 
the present work, we employ the fourth-order scheme, 
with which the convergence is achieved much faster. 
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FIG. 1: Snapshots of the density contour curves for p and the three velocity for «'(= dx^ /dt) in the equatorial plane for model 
A. The solid contour curves are drawn for p — 2i x W^^ g/cm^ (j=l-4) and for 10^^~' g/cm^ (i = 1 ~ 5). The maximum 
density at f = is « 8.86 x 10^^ g/cm^. The blue, cyan, magenta, and green curves denote 10^*, 10^'^, 10^^, and 10^^ g/cm'^, 
respectively. The vectors shows the velocity field {v^,v^), and the scale is shown in the upper right-hand corner. The thick 
(red) circles are the apparent horizons. The 2nd and 3rd panels denote the states after one and one-and-half orbits, respectively. 
The last panel plots the density contours and velocity vectors in the a; = plane at the same time as that for the 8th panel. 
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B. Atmosphere 

Because any conservation scheme of hydrodynamics is 
unable to evolve a vacuum, we have to introduce an arti- 
ficial atmosphere outside the NS. However, if the density 
of the atmosphere was too large, it might affect the or- 
bital motions of the BH and NS and moreover we might 
overestimate the total amount of the rest mass of the 
formed torus surrounding a BH. Thus, we initially assign 
a small rest-mass density as follows: 



Pat r < x{No) 

p^^^l-r/x(No) r>x{No). 



(12) 



Here, we choose pat — Pmax x 10~^ where pmax is the max- 
imum rest-mass density of the NS. With such choice, the 
total amount of the rest mass of the atmosphere is about 
10"^ of the rest mass of the NS. Thus, the accretion of 
the atmosphere onto the NS and BH plays a negligible 
role for their orbital evolution in the present context. In 
the following, we discuss the rest mass of the torus sur- 
rounding a BH. As we show below, the rest mass is larger 
than O.OIMq which is much larger than the atmosphere 
mass. Hence, the atmosphere also plays a negligible role 
for determining the properties of the torus. 



C. General process of tidal disruption 

All the numerical simulations were performed from 
about one and half orbits before the onset of tidal disrup- 
tion to the time in which the accretion time scale of the 
rest mass of material into the BH is much longer than the 
dynamical time scale. The duration of the simulations is 
400-500M in units of the ADM mass. In this section, 
we first describe the general process of tidal disruption, 
subsequent evolution of the BH, and formation of the 
accretion disk (torus) showing the numerical results for 
model A. 

Figure [T] plots the evolution of the contour curves 
for p and the velocity vectors for the three velocity 
v^{= dx^/dt) in the equatorial plane for model A. The 
location of the apparent horizons is shown together. Due 
to gravitational radiation reaction, the orbital radius de- 
creases and then the NS is elongated gradually (panels 
2 and 3). The tidal disruption of the NS by the BH 
sets in at t ^ 4.3 ms (at about one and half orbits; 
panel 3). Soon after this time, the outer part of the NS 
expands outward due to angular momentum transport 
by the hydrodynamic interaction. This material subse- 
quently forms a torus around the BH. However, the tidal 
disruption sets in at an orbit very close to the ISCO and 
hence the material in the inner part is quickly swallowed 
into the BH (4rd and 5th panels) . By the outward angu- 
lar momentum transport, the material in the outer part 
of the NS forms a one-armed spiral arm (5th and 6th 
panels). The spiral arm then winds around the BH, and 
the material, which does not have angular momentum 



large enough to orbit the BH, shrinks and falls into the 
BH (7th panel). A relatively large fraction of the ma- 
terial falls in particular for t « 9-11 ms (see also the 
upper panel of Fig. [5fa)). In this phase, the inner and 
outer parts of the spiral arm collide and thermal energy 
is generated due to the shock heating. The high-density 
part of the spiral arm then expands to form a compact 
torus with the maximum density ~ 10^^ g/cm^ (8th and 
9th panels; see also Fig. [2]). We stopped the simulation 
at t « 12.5 ms, at which the rest mass of the material 
located outside the apparent horizon is ~ O.O92M0. 

We define the specific thermal energy generated by the 
shocks by 



T-r 



(13) 



where the second term in the right-hand side is the poly- 
tropic term, i.e., the specific internal energy in the ab- 
sence of the shocks. Eth is zero in the adiabatic evolution 
(in the absence of the shocks), and hence, eth uiay be 
regarded as the specific thermal energy generated by the 
shocks. Figure [2] (top panel of each) plots eth along x 
and y axes at t = 7.117, 9.070, and 11.02 ms. The ratio 
of the thermal energy to the total internal energy Sth/s 
and the rest-mass density are shown together. This figure 
shows that the generated thermal energy is much larger 
than the polytropic one for most region except for the 
high-density region of the spiral arms which do not expe- 
rience the shock heating. However, during the evolution, 
the spiral arms winds around the BH and the shocked 
region increases. As a result, ffth/e eventually becomes 
approximately equal to unity for all the regions. Figure 
[2] indicates that a hot torus with the density ~ 10^^- 
10^^ g/cm'^ is the final outcome. This is a favorable prop- 
erty for producing a lar ge a mount of neutrinos which may 
drive SGRBs (e.g., HO,!!!, El [H ) . 

Although the EOS adopted in this paper is idealized 
one, we approximately estimate the temperature of the 
torus. Assuming that the specific thermal energy is com- 
posed of those of the gas, photons, and relativistic elec- 
trons (and positrons), we have [H, 



3kT UarT'^ 



£th 



2m 



(14) 



where k, a^, T, and m are the Boltzmann constant, 
radiation-density constant, temperature, and mass of the 
main component of nucleons. For simplicity, we as- 
sume that the gas is composed of the neutron and set 
m = 1.66 X 10"^** g. Because the temperature is higher 
than 6 X 10^ K (see below), we assume that the elec- 
trons and positrons are relativistic to describe Eq. p4p . 
If the density is high, the neutrinos are optically thick 
and contribute to the thermal energy. Although this ef- 
fect is important for p > 10^^ g/cm^ [13, SH, 0, El] , the 
density of the outcome is at most 10^^ g/cm"^, and hence, 
we ignore it. Then, the radiation energy is smaller than 
the gas energy for p > 10"'^^ g/cm'^ and eth/c^ ^ 0.01. 
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FIG. 2: Profiles of eth, Eth/e, and p along x-axis (solid curves) and y-axis (dashed curves) a.t t = 7.117 ms (left), 9.070 ins 
(middle), and 11.02 (right). 




From Eq. (|14p . the temperature is approximately given 
by 

T«.7.2xl0i°(^|^) K. (15) 

Thus, the estimated temperature is between 10^" and 
~ 10"'^'^ K for the torus. (Note that e ~ 0.1 inside the 
apparent horizon, but we do not pay attention to the 
inside of the BH.) 

For clarification of the dependence of the density and 
the thermal energy on the radius of the tidally disrupted 
NSs, i?NS, we show the same plots as Fig. [2]for models B 
and C in Fig. [31 We also plot the density contour curves 
and velocity vectors at t = 11.02 ms for these models in 
Fig. m Figure [3] shows that the value of eth increases 
above 0.001 dominating the specific thermal energy irre- 
spective of models. Figure S] shows that irrespective of 
the models, the outcome after the tidal disruption is a 
compact torus surrounding the BH. Thus, the hot torus 
is the universal outcome irrespective of i?NS as long as it 
is in the range between 12 and 15 km. 



By contrast, the rest- mass density depends sensitively 
on the torus mass (see Table III and next section for the 
final values of the rest mass of the material located out- 
side the apparent horizon which is approximately equal 
to the torus mass). For model B, the torus mass is much 
smaller than that for model A. This is simply because 
the NS radius is smaller and hence the tidal disruption 
sets in at an orbit very close to the ISCO. For model 
C, by contrast, the torus mass is much larger than those 
for models A and B because of the larger NS radius. 
The averaged rest-mass density has a clear correlation 
with the torus mass. For model B, the density is smaller 
than 10^^ g/cm'^ for most region of the torus, whereas 
for model C, the density is larger than 10^^ g/cm^ for 
most region. For such high-density case, the neutrinos 
are likely to be optically thick [4(]| . and the neutrino- 
dominated accretion flow will be subsequently formed. 
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FIG. 4; The same as Fig. 1 but for models B and C ai t — 11.02 ms. The upper two panels plot the contour plots and velocity 
vectors for x-y and y-z planes for model B. The lower two panels show the plots for model C. 



D. Torus mass 

The upper panels of Fig. [5Ua)-(c) plot the evolution 
of the fraction of the rest mass of the material located 
outside the apparent horizon M,.>rAH i^) for models A, 
D, and F, (b) for models B and E, and (c) for models 
A-C. Here, Mr>rAH is defined by 

Mr>r^H = / pau*y/^(fx, (16) 

•^»">rAH 

and u* is the time component of the four-velocity. We 
find that a large fraction of the material is swallowed 
into the BH soon after the onset of merger irrespective of 
the models (see also Table III). However, the final value 
of Mr>r-AH depends strongly on the NS radius; e.g., for 
i?NS = 14.7 km and q w 0.33 (model C), « 12% of the 
material is located outside the BH at the end of the sim- 
ulation, whereas for i?NS = 12-0 km and q w 0.33 (model 
B), only w 2% of the material is located outside the BH 
at the end of the simulation (see Fig. Efc)). The final 
value of Mr>rAH depends also on the mass ratio q. Figure 
[SJa) illustrates that the final value of Mj-^^ah increases 
by ~ 75% for the increase of the value of q from 0.28 to 
0.39. Figure \^h) also illustrates that the final value of 



Mr>r-AH increases by 50% for the small increase of q 
from 0.33 to 0.39. 



Figure in^b) and (c) (in particular (c)) clarifies the de- 
pendence of the infall process of the tidally disrupted ma- 
terial into the BH on the NS radius. With the decrease 
of the NS radii, the accretion rate of the rest mass into 
the BH at the onset of tidal disruption increases steeply. 
Furthermore, the accretion time scale shortens; e.g., for 
model B, the duration of the quick accretion is shorter 
than 1 ms, whereas for model C, it is '-^ 1.5 ms. 



Figure (Hd) plots the values of Afr>rAH the end of 
the simulations as a function of i?NS for given values of 
the mass ratio q. It is clearly seen that the fraction of the 
material located outside the BH steeply decreases with 
the decrease of the NS radius. Extrapolating the results 
for the smaller values of i?NS, we expect that Mr>rAH 
would become close to zero for i?NS ^ H km both for 
q = 0.33 and 0.39. For such case, a very small value of 
the BH mass as Mbh < 3M0 would be required for the 
formation of torus of mass ^ O.OIM0. 
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FIG. 5: The evolution of the rest mass of the material located outside the apparent horizon (upper panel of each figure) and 
the evolution of the area of the apparent horizon in units of IBttM^ (lower panel of each figure) (a) for models A, D, and F 
(for i?NS ~ 13.2 km), (b) for models B and E (for _Rns = 12 km), and (c) for models A, B, and C (for q « 0.33). fdisr in the 
horizontal axis denotes the approximate time at which the tidal disruption sets in, and we determine it to be 4.5 ms for model 
A, 4.2 ms for mode B, 4.3 ms for mode C, 4.2 ms for model D, 3.9 ms for model E, and 5.6 ms for model F, respectively, (d) 
Mr>rAH ^^i*^ '^^ ^he simulatious as a function of -Rns- 



E. Black hole mass and spin 



The lower panels of Fig. [ni^a)-(c) plot the evolution of 
the area of the apparent horizon in units of IGttAI^. The 
area of the BHs is constant before the onset of tidal dis- 
ruption. The value of the area fluctuates only by ~ 0.03% 
in such phase. After the onset of tidal disruption, the 
area quickly increases as a result of the mass accretion, 
and finally, approaches an approximate constant. From 
the final values of the area of the apparent horizon to- 
gether with the estimated BH mass at the final stage, 
MBH.t, the nondimensional spin parameter of the formed 



BHs, a, is approximately derived from 
A 



A: 



1 + VT 



167rM|jj f 



(17) 



To approximately estimate Men.f , we use the following 
equation as in [15| 



Ah 



BH,f 



M -M, 



r>rAH 



Egw, 



(18) 



where ^Eqw is the total radiated energy by gravitational 
waves [52[. The results for Sqw, -^r>rAHJ -^BH.f, and A 
are listed in Table III. The gravitational wave energy is 
primarily carried by the I = \m\ = 2 modes. We found 
that I = ItoI =3 and / = Iml = 4 modes are subdominant 
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TABLE III: The total radiated energy in units of the initial ADM mass M, the total radiated angular momentum in units 
of the initial angular momentum J, the kick velocity, the rest mass of the material located outside the apparent horizon, the 
estimated BH mass, the area of the apparent horizon in units of 167rM|jj f, the ratio of the polar circumferential radius to the 
equatorial one of the apparent horizon (Cp/Ce), and the estimated spin parameters of the final state of the BH. an, af2, and 
flfa are computed from the area of the apparent horizon, the estimated angular momentum and mass of the BH, and Cp/Ce, 
respectively. All the values presented here are measured for the states obtained at the end of the simulations. We note that 
Cp/Ce varies with time by ~ 0.05%. The energy, angular momentum, and linear momentum fluxes of gravitational waves are 
evaluated changing the location of the extraction, and we found that the energy and angular momentum fluxes converge within 
~ 1% error, whereas the error of the linear momentum flux is ~ 10%. 



Model AE/M 


AJ/J 


Vkick (km 




MBH,t(M0) 


Aah 




an 


af2 


af3 


A 0.68% 


11.3% 


19 


0.092 


5.099 


0.9175 


0.9391 


0.550 


0.560 


0.545 


B 0.96% 


12.9% 


13 


0.027 


5.116 


0.9110 


0.9350 


0.570 


0.577 


0.561 


C 0.45% 


8.7% 


18 


0.159 


5.078 


0.9241 


0.9437 


0.530 


0.543 


0.526 


D 0.61% 


10.3% 


17 


0.106 


4.442 


0.8994 


0.9260 


0.597 


0.612 


0.595 


E 0.93% 


13.0% 


19 


0.046 


4.458 


0.8953 


0.9220 


0.612 


0.622 


0.608 


F 0.77% 


12.5% 


4 


0.060 


5.773 


0.9288 


0.9477 


0.514 


0.525 


0.508 



modes; for q w 0.33 (0.39), the total emitted energies 
carried by I — \m\ — 3 and I — \m\ = 4 modes are 5- 
6% (4%) and ~ 1-2% (1%) of that of Z = |to| = 2 modes, 
respectively. 

The nondimensional spin parameter determined from 
A is listed in Table III as an . We find that an is between 
~ 0.5 and ~ 0.6 for all the models, implying that rotat- 
ing BHs of moderate magnitude of the spin arc the final 
outcomes. For the larger mass ratio q, the spin parame- 
ter is larger for the same NS radius. The reason for this 
is that for the larger mass ratio, the initial total angular 
momentum of the system J is larger (see Table I). For 
the smaller NS radius, the spin parameter is larger for 
the same mass ratio. This is simply because the larger 
fraction of the material falls into the BH resulting in the 
spin up. 

To check the validity of the above estimate, we evaluate 
the nondimensional spin parameter employing other two 
methods by which the BH spin may be approximately 
calculated. In the first method, the spin is evaluated 
directly from the angular momentum of the BH defined 

by 

JbHJ = J — Jr>rAH ^ JgW, (19) 

where Jqw is the total radiated angular momentum by 
gravitational waves and Jr>rAH is th^ angular momen- 
tum of the material located outside the apparent horizon, 
which is defined by 

Jr>rAH = / pahu*u^py/j(fix. (20) 

Here, is the ip component of the four-velocity, and h is 
the specific enthalpy defined by 1+e+P/ p. JsH.f exactly 
gives the angular momentum of the material in the ax- 
isymmetric and stationary spacetime. In the late phase 
of the merger, the spacetime relaxes to a quasistationary 
and nearly axisymmetric state. Thus, we may expect 
that Jr>rAH '^iil provide an approximate magnitude of 
the angular momentum of the torus. 



From JBH,f and MeH.f, we define the nondimensional 
spin parameter by af2 = JBH,f/AfBH f (^^^ Table HI). It is 
found that a{2 is systematically larger than an but these 
two values, independently determined, still agree within 
w 2.5% disagreement. This suggests that both quantities 
denote the BH spin within the error of ~ 2.5% (^ 0.015). 
We note that because the relation Ofi < at2 holds system- 
atically, this disagreement would primarily result from 
the systematic error associated with the approximate def- 
initions of the spin. For example, we use the rest mass 
Mr>rAii to estimate the disk mass in Eq. (llSp . ignoring 
the gravitational binding energy between the BH and the 
torus. Also, Eq. pp)) may include an error because the 
torus is neither exactly axisymmetric nor stationary. 

In the second method, we use the ratio of the polar 
circumferential length to the equatorial one of the BH, 
Cp/Ce- We compare this ratio to that of the Kerr BH, 
and estimate the spin parameter. The results are listed in 
Table HI as ots. Again we find that this agrees with other 
two with a good accuracy (within ~ 1.5% disagreement 
with an ) ■ We note that the axial ratio varies with time by 
~ 0.05%. This implies that the error for the estimation 
of at3 is ~ 0.5% for ata = 0.5-0.6. We also note that the 
relation between Cp/Ce and the spin parameter in the 
presence of the torus is slightly different from that of the 
Kerr BH [ij . However, the systematic error is negligible 
because the torus mass is only % 3% of the BH mass. 

From these results, we conclude that the spin param- 
eter obtained by three different methods agree within 
~ 2.5% error. We note that all three methods can only 
approximately determine the spin, and the systematic 
error associated with their definition should be included. 
Actually, the systematic relation af2 > an > ots holds. 
Hence, the error of ~ 2.5% results not only from the 
numerical error but also from the systematic error. 

The spin parameter of the formed BHs is smaller than 
the initial value of J/A'P of the system by 14-20%. The 
primary reason for this in the case that most of the NS 
matter falls into the BH, i.e., for i?NS = 12 km, is that 
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gravitational waves carry away the angular momentum 
by ~ 13% of J (see Table II). By contrast, for the case 
that AIj'^j-j^^ is large, the angular momentum of the torus 
occupies a large fraction of J; e.g., for model C, in which 
Mr>r-AH ~ 0.16Mq, Jr>rAH about 15% of J. In such 
case, the angular momentum of the torus is about twice 
larger than that carried by gravitational waves. 



F. Gravitational waves 

Figure [U^a) and (b) plot gravitational waveforms ob- 
served along the z-axis as a function of the retarded time 
for model A. The retarded time is approximately defined 
by 



irct = t- D - 2M\n{D/M), 



(21) 



where D is the distance to the observer. The upper panel 
of Fig. mja) plots + and x modes of gravitational waves 
and the lower one the Newman-Penrose quantity. The 
gravitational waveforms are obtained from twice time in- 
tegration of the Newman-Penrose quantity. Because the 
I = \m\ =2 mode is dominant in the observation along 
the z-axis, we determine the waveforms only from this 
mode. 

From the values of and h-^ , the amplitude of grav- 
itational waves at a distance D is evaluated 



/igw « 5.0 X 10" 



h\ 



0.2 

100 Mpc 



D 



)( 



Afn 



5.2Af, 



(22) 



Figure EJa) implies that the maximum amplitude at a 
distance oi D = 100 Mpc is « 4 x 10"^^ for model A. 

For trot ^ 4.5 ms, the inspiral waveforms are seen: 
The amplitude increases and characteristic wavelength 
decreases with time. The wavelength at the last inspi- 
ral orbit is 0.8 ms, indicating that the orbital period 
at the complete tidal disruption is ~ 1.6 ms, i.e., the 
angular velocity is ~ O.IOA/"^. This value coincides ap- 
proximately with that of the ISCO, but does not agree 
with Vt calculated from Eq. ([4]) (see Sec. I for discus- 
sion) . This indicates that the tidal disruption may set in 
approximately at the orbit predicted by the study of the 
quasicircular orbits but completes near the ISCO. 

For 4.5 ms ^ irct ^ 5.5 ms, the amplitude of gravita- 
tional waves decreases quickly. In this phase, the tidal 
disruption and resultant quick accretion of the material 
into the BH proceed (see Fig. [5|). Thus, gravitational 
waves could be emitted both by the matter motion and 
by the quasinormal mode (QNM) oscillation of the BH. 
The characteristic oscillation period in this phase is ~ 2.9 
kHz. This value is slightly smaller than the value pre- 
dicted from the perturbation study for the fundamental 
QNM of BH mass MsH^f ~ 5.10Mo and spin a = 0.55. 



Here, the perturbation study predicts the frequency and 
damping time scale as (Zsj 



td 



3.23 



A^BH.f 



, lOAfQ 
2(1 -a)^"-45 



[1 - 0.63(1 - af -"] kH^33) 



(24) 



which gives /qnm ~ 3.19 kHz. This indicates that gravita- 
tional waves are primarily emitted by the matter motion 
in this phase; the damping of the amplitude is not due to 
the QNM damping but to the fact that the compactness 
and degree of nonaxisymmctry of the matter distribution 
decrease by the tidal disruption. 

For t > 5.5 ms, the rapid accretion stops and the 
BH approaches a nearly stationary state (see Fig. [5ja)). 
Hence, for t-^ct ^ 5.5 ms, the ring-down waveforms of the 
QNM should be seen. Because gravitational waves are 
also emitted by the matter moving around the BH, the 
waveforms are modulated and the waveforms associated 
only by the QNM are not clearly seen. However, '04 with 
5.4 ms < tj-ct ^ 6 ms can be fitted by the damping wave- 
forms of the fundamental QNM fairly well. Figure EJb) 
plots the I = \m\ = 2 mode of 1/^4 and a fitted waveform 
as 



Ae-'/'^ sin(2^/q„„,i + (5). 



(25) 



where A and 5 are constants. For the fitting, we choose 
the BH mass and spin as 5.10A/q and a = 0.55. In this 
case, /qnm = 3.19 kHz and td — 0.286 ms, respectively 
[4^ . This figure shows that gravitational waves are pri- 
marily characterized by the fundamental QNM, and the 
gravitational wave amplitude of the QNM is much smaller 
than that at the last inspiral orbit. The possible reason 
for this small amplitude is that the QNM is excited only 
incoherently by the infall of the tidally disrupted non- 
compact material into the BH. 

The characteristic feature of the gravitational wave- 
forms emitted after the tidal disruption sets in is that 
the amplitude damps quickly even in the formation of 
the nonstationary torus of mass ^ 0.1A/q. This is due to 
the fact that the compactness and degree of nonaxisym- 
metry of the torus decrease in a very short time scale 
{'-^ 2-3 ms). Moreover, the gravitational wave amplitude 
of the QNM is very small. This implies that the am- 
plitude of the Fourier spectrum steeply decreases in the 
high-frequency region for / > 1/0.8 ms « 1.2 kHz (see 
Sec. UVUI cf. Fig. [T]). 

Gravitational waveforms are qualitatively similar for 
other models. However, quantitative differences due to 
the difference of the NS radius are clearly seen. Figure [H] 
(c) plots gravitational waves for model C (solid curves) 
together with those for model A (dashed curves). For 
t-cct ^ 4 ms, the inspiral waveforms for two models agree 
very well, because gravitational waves at such phase are 
primarily determined by the total mass and mass ratio. 
However, for ij-et ^ 4 ms, the waveforms disagree be- 
cause the amplitude for model C starts decreasing (com- 
pare the waveforms for 4.5 ms ^ ijct ^ 6 ms). This 
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FIG. 6: (a) upper panel: + and x modes of gravitational waveforms observed along the z-axis for model A (solid and dashed 
curves, respectively), tret denotes the retarded time [see Eq. (I21|l ] and Mo is the total mass defined in the caption of Table 
I. The amplitude at a distance of observer can be found from Eq. (|22|) . (a) lower panel: The real and imaginary parts of tp4, 
(solid and dashed curves), (b) The real (upper panel) and imaginary (lower panel) parts of tpi for model A together with the 
hypothetical curves of gravitational waves associated with the QNM (dotted curves) for 5 ms ^ tret i$ 6.5 ms. (c) h+ (upper 
panel) and Re(-i/)4) (lower panel) for models A (dashed curve) and C (solid curve), (d) The real (solid curve) and imaginary 
(dashed curve) parts of ^pi for model B (upper panel). The waveforms for 4.8 ms ^ tret ^ 6 ms are magnified in the lower panel 
together with the hypothetical curves of gravitational waves associated with the QNM (dotted curves). For all the plots, only 
the I = Iml = 2 modes are taken into account. 



earlier decrease results from the onset of tidal disruption 
of the NS at a larger orbital separation due to the larger 
NS radius for model C (see Eq. [4]). Soon after the on- 
set of tidal disruption, the material spreads around the 
BH at a relatively large orbital separation. Because the 
compactness and degree of nonaxisymmetry decrease as 
the tidally disrupted material spreads, the amplitude for 



model C quickly damps. Furthermore, the fraction of the 
material coherently falling into the BH is not large. As 
a result, the amplitude of the QNM is negligible (com- 
pare ■04 of models A and C; because it is smaller than 
the amplitude of the wave modulation due to numerical 
error or matter motion, we cannot extract the QNM for 
model C). 
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By contrast, the QNM is clearly seen for the case 
that the NS radius is smaller. Figure [DJd) plots the 
I = \m\ — 2 mode of ■04 for model B. In this case, the 
tidal disruption sets in at an orbit closer to the ISCO 
than that for models A and C. Consequently, a larger 
fraction of the material coherently falls into the BH, and 
relatively coherently excites the QNM. In the lower panel 
of Fig. [HJd), we again fit the waveforms by Eq. (pS)) with 
/qnm = 3.23 kHz and td = 0.289 ms assuming that the 
BH mass and spin are 5.12M0 and 0.57, respectively. 
Figure [6] (d) shows that the waveforms for ij-et ^ 5 ms 
agree well with the hypothetical fitting formula. A dis- 
tinguishable feature is that the amplitude of the QNM is 
by a factor of ~ 2 larger than that for model A. However, 
this amplitude is still much smaller than the amplitude 
at the last inspiral orbit. 

The kick velocity is also estimated from the total linear 
momentum carried by gravitational waves. We find that 
it is < 20 km/s for aU the models (see Table III) [53l |. 
This value is much smaller (by one order of magnitude) 
than that in the case of the BH-BH merger [38|. The 
reason for this is that the gravitational wave amplitude 
damps quickly during the merger due to the tidal dis- 
ruption, and the amplitude associated with the QNM is 
very small. Actually, a PN study [i^ illustrates that the 
kick velocity is « 20 km/s for q = 1/3 if the kick only by 
gravitational waves from the inspiral motion is taken into 
account. This value agrees fairly with our results. In the 
BH-BH merger, gravitational waves associated with the 
QNM have a large amplitude and as a result, the linear 
momentum flux is significantly enhanced. If the NS es- 
capes the tidal disruption during the merger (i.e., for the 
case that the BH mass is large enough or the NS radius 
is small enough), a large kick velocity of order 100 km/s 
might be induced. 



G. Gravitational wave spectrum 

To determine the characteristic frequency of gravi- 
tational waves emitted during the tidal disruption, we 
compute the Fourier spectrum of gravitational waves of 
Z = |r7i| = 2 modes defined by 



D l\h+if)\^ + \h.{fW 



Mo 



where 



h+{f) 



(26) 

(27) 
(28) 



In the present simulations, we followed the inspiral phase 
only for about one and half orbits before the onset of 
tidal disruption. Thus, the Fourier spectrum for the 
low-frequency region is absent if we perform the Fourier 
transformation for the purely numerical data (see the 



dashed curve of Fig. [HJa)). To artificially compensate 
the Fourier spectrum for the low-frequency region, we 
combine the hypothetical waveforms computed from the 
third PN approximation for two structureless, inspiraling 
compact objects |6(|. Specifically we determine the evolu- 
tion of the orbital angular velocity Q. using the so-called 
Taylor T4 formula (e.g. [131 for a detailed review) by 
solving 
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where x = [MQD.{t)Y/^ is a function of time [5j], 77 is 
the ratio of the reduced mass to the total mass Mq, and 
7e = 0.577 • • • is the Euler constant. Then, gravitational 
waveforms are determined from 

h+{t) = ^^^^^(a;) cos[$(0 H- 5], (30) 
/ix(t) = ^^A(x)sin[$(t)-HJ], (31) 

where A{x) is a nondimensional function of x for which 
A{x) ^ 1 for a; ^ 0, (5 is an arbitrary phase, and 



$(t) = 2 j n{t)dt. 



(32) 



The 



For A{x), we choose the 2.5 PN formula (e.g., [33]). 
PN waveforms are calculated for Mofi > 0.005. 

In this method, it is not clear at which frequency we 
should combine the PN waveforms with the numerical 
waveforms. In the present paper, we match the two wave- 
forms at Mon « 0.04 (at / « 517(A'fo/5M0)-i Hz). The 
Fourier spectrum for the late inspiral phase depends on 
this matching frequency, in particular, around the match- 
ing frequency ~ 500-700 Hz. However, the spectrum 
around the frequency for / > /tidal, for which the tidal 
disruption proceeds, depends very weakly on the match- 
ing frequency. Thus, the present method is acceptable 
for studying the Fourier spectrum during the tidal dis- 
ruption. 

Figure[7Ka) plots h{f)f for model A. We also plot the 
spectra for the analytical third PN waveforms (dotted 
curve) and for the purely numerical waveforms with no 
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FIG. 7: (a) Spectrum of gravitational waves h{f)f for model A (solid curve). The dashed and dotted curves are the spectra 
for the purely numerical (short-term) data and for the analytic third PN waveform, (b) the same as (a) but for models A, 
D, and F. The vertical dotted lines denote the expected frequencies at which the tidal disruption sets in for models D and F; 
/tidal ~ 0.84 kHz for model F and 0.88 kHz for model D. (c) the same as (b) but for models D and E. The vertical dotted line 
denotes the expected frequency at which the tidal disruption sets in, /tidal ~ 1-0 kHz. (d) the same as (b) but for models A-C. 
For (b)-(d), the solid circles denote the location of /cut. 



matching (dashed curve). For nM^f <C 1, the spectra of 
the matched and PN waveforms behave Q 

Kf)f-\[^n"'{^M,f)-'/'. (33) 

We note that the effective amphtude defined by the av- 
erage over the source direction and the direction of the 
binary orbital plane is 

_ 2 fe(/)/Afo 

"off = 

V 0.1 )\bMQ) 

(34) 

V 100 Mpc J ^ ' 



Thus, the effective amplitude is « 10 at / w 1 kHz 
for D = 100 Mpc and Mq = 5Mq. 

It is found that the spectrum for the purely numeri- 
cal data agrees with the spectrum of the matched data 
for / > 800 Hz besides the spurious modulation for the 
spectrum of the purely numerical data which results from 
the incomplete data sets for the inspiral phase. Here, 
/ = 800 Hz is slightly smaller than /tidal predicted from 
Eq. (O (see also Table IV). This illustrates that the 
matching does not affect the global shape of the spec- 
trum for / > /tidal- 

We still see the modulation and dip of the spectrum for 
the matched data in the frequency band between ~ 500 
and 800 Hz. The possible physical reason for the dip at 
/ 1 kHz is that the radial-approach velocity of the BH 
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and NS steeply increases near the ISCO (/ ~ 1.2 kHz), 
and hence, the integration time for gravitational waves 
decreases. As a result, the effective amplitude at such 
high-frequency band slightly decreases. This feature is 
well-known for the NS-NS merger [47[. However, there 
are also the possible unphysical reasons as follows: (i) the 
orbit of the BH-NS has a nonzero eccentricity because the 
initial data is not exactly the circular orbit [48] , and (ii) 
the matching of the numerical and PN data induces a 
systematic error. Thus, the spectrum for this band may 
not be very accurate. We do not pay attention to this 
low-frequency band in the following. 

The noteworthy feature of the spectrum is that the 
spectrum amplitude does not damp even for /tidal < / 
1.3/tidai = /cut (see the sohd circles in Fig. WUo)-{d) 
for the approximate location of /cut)- The reason for 
this is that the NS orbiting near the ISCO has a radial- 
approach velocity of order ~ 10% of the orbital velocity 
due to the radiation reaction of gravitational waves, and 
hence, the NS is not immediately tidally disrupted at the 
orbit of / = /tidal, although the tidal disruption of the 
NS hkely sets in at / = /tidal- As a result, the inspi- 
ral orbit is maintained for a while even inside the pre- 
dicted tidal disruption orbit and gravitational waves of a 
large amplitude are emitted for / > /tidal- Because the 
tidal disruption completes at / « /cut > /tidal, it is not 
straightforward to determine /tidal from the spectrum of 
gravitational waves. 

We note that our finding of /cut > /tidal is the unique 
feature for the BH-NS binary. If the companion of the BH 
is not as compact as the NS (e.g., for the BH- white dwarfs 
binary), /cut would be approximately equal to /tidal be- 
cause the radial-approach velocity due to the gravita- 
tional radiation reaction is negligible and hence the tidal 
disruption would complete at the predicted tidal disrup- 
tion orbit. 

Because the amplitude of gravitational waves quickly 
decreases after the tidal disruption completes, the spec- 
trum also sharply falls for the high-frequency band; i.e., 
above a cut-off frequency /cut, the spectrum amplitude 
steeply decreases. Figure ^a) shows that for model 
A, the value of /cut is approximately 1.16 kHz (i.e., 
Mo^ « 0.094). This confirms the finding in Sec. HVFl 
about the quick decrease of the gravitational wave am- 
plitude occurs near the ISCO in this case. 

To see the dependence of the spectrum shape on the 
mass ratio and NS radius, we plot Fig. [Hb)-(d). Figure 
[7]Jb) and (c) compare the spectra for models A, D, and 
F (i?NS = 13.2 km) and for models B and E (i?NS ~ 
12.0 km), for which the NS radii are identical each other, 
whereas the mass ratio q is different. It is found that the 
spectrum shape depends weakly on the mass ratio. The 
reason for this is that the gravitational wave frequency at 
which the tidal disruption occurs is determined primarily 
by the radius and mass of the NS. For models A, D, and 
F, /cut ~ 1.1-1.2 kHz, and for models B and E, /cut ~ 
1.4 kHz (see Table IV). Obviously, for the smaller radii, 
the value of /cut is larger because the tidal disruption 



TABLE IV: The expected frequency of gravitational waves at 
the tidal disruption /tidai, the peak frequency of gravitational 
wave spectrum near the sharp cut-off in the high-frequency 
region /cut, and the ratio /cut//tidai. 

Model /tidal (kHz) /cut (kHz) /cut //tidal 



A 


0.856 


1.16 


1.36 


B 


0.997 


1.41 


1.41 


C 


0.736 


0.92 


1.25 


D 


0.877 


1.14 


1.30 


E 


1.021 


1.40 


1.37 


F 


0.840 


1.09 


1.30 



completes for the smaller orbital separation. 

We also calculate the ratio of /cut to /tidal (see Table 
IV). For models A, D, and F, the ratio is ~ 1.3, whereas 
for models B and E, it is slightly larger, ^ 1.4. The reason 
is that the compact NS has the larger radial-approach 
velocity at / = /tidal, and hence, the tidal disruption 
completes at a smaller orbital separation. 

Figure (Tljd) compares the spectra for models A-C, 
for which the mass ratio q is approximately identical, 
whereas the NS radii are different. We find that the 
value of /cut depends strongly of the NS radius and is 
smaller for larger NS radii, indicating that the tidal dis- 
ruption sets in at a larger orbital separation. The ratio 
/cut //tidal is also Smaller for larger NS radii (see Table 
IV). The reason for this is that for the larger NS, the 
tidal disruption sets in at a larger orbital separation, and 
hence, the radial-approach velocity is smaller. However, 
even for model C, the ratio is « 1.25. This implies that 
gravitational waves of a high amplitude are emitted well 
inside the orbit of the onset of tidal disruption, even for 
a less compact NS of nearly upper-limit radius 15 km. 

Finally, we note the following: Figure [7] and Equation 
([M)) show that the effective amplitude at / = /cut ^ 1 
kHz is « 10^^^ for the typical distance and total mass as 
D — 100 Mpc and Mq — 5Mq. Even for the optimistic 
direction of the source and its binary orbital plane, the 
effective amplitude is at most « 2.5 x 10~^^. The de- 
signed sensitivity of the advanced LIGO is ~ 3 x 10~^^ 
[4g |. This implies that it will not be possible to detect 
gravitational waves during the tidal disruption by the de- 
tectors of standard design. To detect gravitational waves 
at such high frequency, the detectors of special instru- 
ment (e.g., resonant-side band extraction [50]) which is 
sensitive to high-frequency gravitational waves is neces- 
sary. 



V. SUMMARY 

We have presented the numerical results of fully gen- 
eral relativistic simulation for the merger of BH-NS bina- 
ries, focusing on the case that the NS is tidally disrupted 
by a nonspinning low-mass BH of A/bh ~ 3.3-4.6M0. 
The F-law EOS with F = 2 and irrotational velocity field 
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are employed for modeling the NSs. To see the depen- 
dence of the tidal disruption event on the NS radius, we 
choose it in a wide range as i?NS = 12.0-1.47 km whereas 
the NS mass is fixed to be « l.SM©. The resulting mass 
ratio q = Mns/A/bh is in the range between 0.28 and 

0. 4. 

As predicted by the study for the quasicircular states 
0) for all the chosen models, the NS is tidally dis- 
rupted at the orbits close to the ISCO. The BH of the 
spin of 0.5-0.6 is formed and a large fraction of the ma- 
terial is quickly swallowed into the BH, whereas 2-12% 
of the material forms a hot and compact torus around 
the BH. The resultant mass and density of the torus de- 
pend strongly on the mass ratio q and NS radius i?NS, 
in particular on i?NS- For i?NS = 12-0 km, the torus 
mass is at most O.OSM© even for the large value of the 
mass ratio q ~ 0.39. For i?NS = 14.7 km, by contrast, 
the torus mass is w O.IGMq for q « 0.33. Extrapolat- 
ing the results in this paper, the torus mass would be 
smaller than O.OIM© for i?NS = 11 km even for q ^ 0.4. 
The stiff nuclear EOSs predict the radius of w 11-12 
km for A/ns = 1.3~1.4Mo [29|]. This suggests that the 
torus mass is likely to be ^ O.IMq even for the BH mass 
3-4M0. This point should be more rigorously clarified 
for the future simulation employing the realistic nuclear 
EOSs. 

For the optimistic cases in which i?NS ^ 13 km or 
of g > 0.4, the torus mass can be > 0.05Mq. As 
we found from the value of the eth, the resultant torus 
likely has high temperature as lO-'^'^-lO"'^"'^ K. This suggests 
that such outcomes are promising candidates for driving 
SGRBs. According to the latest simulations for the hot 
and compact torus around the BH by Setiawan et al. [iH , 
the total energy deposited by neutrino-antineutrino an- 
nihilation is ~ 1049(Mtorus/0.01MQ) ergs for O.OIMq < 
Aftorus ^ O.IMq and a = 0.6. Here, Mtorus is the initial 
torus mass for the simulation. Assuming that the con- 
version rate to the gamma-ray energy is 10% [5l|, the 
total energy of SGRBs is ~ 10'^*(Mtorus/0.01Mo) ergs 
according to their numerical results. This indicates that 
if the NS radius is > 13 km with Mbh = 3.3-4Mo, the 
SGRB of the total energy ~ lO''^ ergs may be driven. 
By contrast, for /Jns ^ 12 km with A^bh ^ 3.3A'f0, the 
total energy is likely to be at most several xlO'** ergs; 

1. e., only the weak SGRBs can be explained. We note 
that the final spin parameter depends on the initial spin 
of the BH. In the merger of the rapidly rotating BH and 
NS, the resuling spin of the BH will be close to unity. For 
such BHs, the total deposited energy will be enhanced by 
the spin effects (43j . 

Gravitational waves are also computed. It is found 
that the amplitude of gravitational waves damps during 
the tidal disruption, and when the tidal disruption com- 
pletes, the amplitude becomes much smaller than that at 
the last inspiral orbit. Although we find that the wave- 
forms in the final phase of the tidal disruption are char- 
acterized by the QNM ringing, its amplitude is much 
smaller than that of gravitational waves at the last in- 



spiral phase, in particular, for the large value of i?NS- 
The reason for this is that the NS is tidally disrupted 
before plunging into the BH, and hence, a significant ex- 
citation of the QNM by the coherently infalling material 
is not achieved in contrast to the BH-BH merger (e.g., 
[H [H, HO, El). For the case that the NS is compact 
with i?NS = 12 km, the tidal disruption occurs at an 
orbit very close to the ISCO and a large amount of the 
material falls fairly coherently into the BH, resulting in a 
relatively high amplitude of the QNM. However, the am- 
plitude is still much smaller than that at the last inspiral 
orbit. 

The Fourier spectrum of gravitational waves is also 
analyzed. Because the amplitude of gravitational waves 
quickly damps during the tidal disruption, the spectrum 
amplitude steeply decreases above a cut-off frequency 
/cut- The noteworthy feature is that the cut-off frequency 
does not agree with the predicted frequency at which the 
tidal disruption sets in, /tidal- The reason for this is that 
the NS in the compact binaries has a radial-approach ve- 
locity of order ~ 10% of the orbital velocity and is not 
immediately tidally disrupted a,t f — /tidal- We have 
found that /cut //tidal is in the range between 1.25 and 
1.4 in our results depending primarily on the NS radius 
-Rns; this value is larger for the smaller values of ii!NS- 
These results imply that it is not straightforward to de- 
termine /tidal from gravitational wave observation. If the 
dependence of /cut //tidal on the NS radius and mass ra- 
tio is clarified in detail by the numerical simulation, /tidal 
may be inferred from /cut- If such a relation is found, /cut 
will be useful for determining the properties of NSs such 
as the radius and density profile. For this purpose, a de- 
tailed simulation taking into account the nuclear EOSs 
is necessary. We plan to perform such simulation in the 
next step. 
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APPENDIX A: COMPARISON WITH THE 
LOW-RESOLUTION RESULTS 

In this appendix, we present the results for models 
Al, A2 and C2 for which the physical parameters of the 
initial conditions are the same as those for models A and 
C but the simulations were performed with the poorer 
grid resolutions (see Table II). 
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FIG. 8: The same as Fig. [^c) but for models A, C, Al, A2, 
and C2. We choose the values of fdisr = 4.30 ms, 4.17 ms, 
and 4.50 ms for A2, Al, and C2, respectively. For models A 
and C, see the caption of Fig. |5l 



Figure [8] plots the evolution of the rest mass of the 
material located outside the apparent horizon Mryrj^n 
and the area of the apparent horizon in units of IGttM^ 
for models Al, A2, and C2. For comparison, the results 
for models A and C are shown together. 

We find that the results for models Al and A2 agree 
qualitatively well with those for model A. Because of the 
difference of the grid structure for covering the BH and 
NS and/or because of the larger numerical dissipation for 
the low-resolution runs, the time of the onset of tidal dis- 
ruption tdisr slightly disagrees among three models, but 
the difference is only « 0.2-0.3 ms 10% of the orbital 
period at the last inspiral orbit). The value of Mryr ah 
(Aah/ISttM^) is systematically overestimated (under- 
estimated) for the lower-resolution runs. However, the 
magnitude of the disagreement is not large. For model 
A2 (Al), the relative errors of MryrAa ^ah to the 
values for model A is ~ 7% (- 4%) and - 0.3% (~ 0.3%), 
respectively. Furthermore, the values systematically con- 
verge, although the order of the convergence is not ex- 
actly specified because the difference among three results 
is small. (Note that the order of the convergence in the 
presence (absence) of shocks should be the first (third) or- 
der for the hydro part whereas that for the gravitational 
field is fourth order, and thus, the order of the conver- 
gence is not simply determined). For determining the 
order of the convergence, it is necessary to perform simu- 
lations with a better accuracy. However, it is not feasible 



we find that the errors for Mr>rAH t 



idisr = 7 ms is 

- O.O19M0 (~ 22%), 0.024Mq (- 30%), and O.O3OM0 
(~ 35%) for models A, Al, and A2, respectively. For 
Aah/ISttA/^, the relative errors at i — idisr = 7 ms are 

- 0.3%, 0.5%, and 0.5% for models A, Al, and A2. 

We find that the results for model C2 agree well with 
those for model C. In this case, the time of the onset 
of tidal disruption tdisr agrees within « 0.04 ms (~ 2% 
of the orbital period at the last inspiral orbit). Further- 
more, the value of Mr>rAH ^'^d the area of the apparent 
horizon agree within ~ 3% and ~ 0.1% errors, respec- 
tively. Assuming the second-order convergence, the ex- 
act values of Mr>rAH ^'^d ^ah/IBttM^ at t — ^disr = 8 
ms are determined, and by comparison with these val- 
ues, we find that the errors of these values for model C 
(C2) are 4% (6%) and 0.2% (0.4%), respectively. For 
model C2, the grid spacing for the inner computational 
domain is Ax = Mp/12, with which the major diameter 
of the NS is covered by the 50 grid points (see Table II). 
Thus, with such setting, the orbits of the BH and NS 
are computed with an accuracy enough at least for the 
qualitative study. 

Figure [H plots and (a) for models A and A2 and 
(b) for models C and C2. This shows that the waveforms 
for the high and low-resolution runs agree qualitatively 
well. For the inspiral phase, the difference of the wave- 
forms for the two different resolutions is very small. In 
particular, for models C and C2, the amplitude at a given 
time agrees each other within ~ 3% error for tj-ot ^4.5 
ms. During the tidal disruption phase, the amplitude is 
systematically smaller for the low-resolution runs. The 
possible reason for this is that the compactness of the 
NS is more quickly lost with the lower resolution, result- 
ing in the less coherent excitation of gravitational waves. 
However, the difference of the amplitude is still at most 
~ 10% for model C2. For model A2, the error is larger in 
particular for 4.5 ms ^ iret ^ 5 ms, although the wave- 
forms agree qualitatively with those for model A. The 
primary source of the error is the phase difference caused 
by the fact that for the lower resolution, the tidal disrup- 
tion sets in earlier. During the tidal disruption, the wave 
amplitude for model A2 is by a factor of ~ 2 smaller than 
that for model A. This indicates that the simulation of 
the poor resolution significantly underestimates the am- 
plitude during the tidal disruption phase. Because of the 
underestimation of the amplitude, the total energy and 
angular momentum emitted by gravitational waves are 
underestimated by ~ 5% for model C2 and by ~ 10% for 
model A2 in comparison with models C and A. 



[1] I. H. Stairs, Science 304, 547 (2004). 

[2] R. Voss and T. M. Taulis, Mon. Not. R. Astro. Soc. 



342, 1169 (2003); K. Belczynski, T. Bulik, and B. 
Rudak, Astrophys. J. 571, 394 (2002); V. Kalogera et 



19 



Q 
+ 

xi 



a 
si 



(a) 




Q 

+ 



a 

X 



(b) 



0.1 


-0.1 

0.1 


-0.1 





2 4 6 

tret (ms) 
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